Crystallization in a dense suspension of self-propelled particles 
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Using Brownian dynamics computer simulations we show that a two-dimensional suspension of 
self-propelled ("active") colloidal particles crystallizes at sufficiently high densities. Compared to 
the equilibrium freezing of passive particles the freezing density is both significantly shifted and 
depends on the structural or dynamical criterion employed. In non-equilibrium the transition is 
accompanied by pronounced structural heterogeneities. This leads to a transition region between 
liquid and solid in which the suspension is globally ordered but unordered liquid-like "bubbles" still 
persist. 
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Recently, the collective dynamics of self-propelled 
("active") particles has become a topic of intense re- 
search [1] [2] resulting in a wealth of new non-equilibrium 
phenomena like swarming [31 0] , clustering ^-(7^ and ac- 
tive swirling [S]. These phenomena have been observed 
both in dense bacterial solutions [9] and in artificial 
microswimmers [TO]. Excellent model systems for self- 
propelled particles are colloidal suspensions, where the 
motility of colloidal particles can be achieved and steered 
by magnetic beads acting as artificial fiagella [11] , by cat- 
alytic reactions at Janus-particles [T2] , or by laser-heated 
metal-capped particles |13 . 

The purpose of the present Letter is to show that self- 
motile interacting colloidal particles in two dimensions 
still freeze into a crystalline lattice displaying long-ranged 
orientational order despite the fact that energy is in- 
jected incessantly. We explore the nature of this non- 
equilibrium transition by Brownian dynamics computer 
simulations of a Yukawa model of self-propelled particles. 
We use a minimal model without explicit alignment of 
particle orientations. In equilibrium, i.e. in the absence 
of self-propagation, freezing and melting of colloidal sus- 
pensions is well understood [14[ . But even for passive 
particles it is known that freezing is seriously affected and 
changed under non-equilibrium conditions, e.g. in a time- 
oscillatory external force field or in shear flow [IB] . 
Recently it has also been shown that active matter can 
reach steady states with frozen fluctuations [17]. 

For self-propelled particles we find that the freezing 
transition is largely shifted relative to its equilibrium lo- 
cation. This shift cannot be explained by a simple scaling 
using the concept of an effective temperature [IB] ; quite 
in contrast to sedimentation profiles of suspensions |19j or 
the long-time diffusion of single propelled particles [20] . 
Rather, the transition points based on different criteria 
for melting and freezing, which agree in equilibrium, di- 
verge. In particular, the dynamical Lindemann-like melt- 
ing [m [22] and freezing criteria [23l [24] define a tran- 
sition region between liquid and solid characterized by 
inhomogeneities of the orientational order parameter. 

We study a suspension of N self-propelled particles 
moving in two dimensions and immersed in a solvent. 



Even though the particles are driven we assume that the 
solvent remains in equilibrium at the well-defined tem- 
perature T. The overdamped motion of the ith particle 
is described through 

r, = -V,;7 + /e,+e.. (1) 

The noise models the stochastic interactions with the 
solvent molecules. It has zero mean and correlations 
{ii{t)i3 {i')) = 2(5.yl(5(t - t'), where 1 is the identity 
matrix. Throughout the paper we employ dimensionless 
quantities and measure energy in units of k^T , length in 
units of and time in units of {pDq)~^ . Here, p is 

the number density and Dq is the bare diffusion coeffi- 
cient. Particles interact pairwise through the repulsive 
Yukawa potential 

uir) = r— (2) 
r 

with screening length A^^ and dimensionless coupling pa- 
rameter r = Voy^/fceT, where Vq is the bare poten- 
tial strength. The total potential energy then becomes 
U — J2i<j "^il^i ~ III addition to the conservative 

force due to J7 a constant force / propels every particle 
in the direction 

^ ( ) ' = 2A.%5(t - t'). (3) 

In the minimal model studied here we assume that these 
particle orientations undergo free diffusion without ex- 
plicit alignment. For spherical particles with diameter a 
the rotational diffusion coefficient is D,- = 3Do/a'^. 

We perform Brownian dynamics simulations for TV = 
1936 particles using periodic boundary conditions. Com- 
mensurable box dimensions L^/Ly = 2/\/3 are chosen 
such that the suspension can crystallize into the hexag- 
onal crystal without any defects. We fix the rotational 
diffusion coefficient to Dr = 3.5 and the inverse screening 
length to A = 3.5; leaving F and / as variable parame- 
ters. The time step for updating particle positions is 
At = 10^^, while particle-particle interactions are cut 
off after an inter-particle distance of 7/A = 2. 
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FIG. 1: (Color online) Cooling (solid lines) and melting curves (dashed lines) for (A) the orientational order parameter i/)6 and 
(B) the long-time diffusion coefficient D vs. the potential strength F for selected driving forces /. The crossings with the dashed 
horizontal lines define the position of the structural transition Fg (i/^e = 0.45) and the dynamical freezing FJ^ {D — 0.086), 
respectively. (C) Phase diagram in the f-F plane. The symbols mark the numerically estimated dynamical freezing line F^ 
and melting line F£ (see main text for definition). The thick dashed line indicates the structural transition Fg. Also plotted 
are the ipe = 0.67 and tpQ = 0.8 "iso-structure" lines along which ipe is constant. 



We simulate cooling and melting runs for forces ^ 
/ ^ 15. For the cooling runs we use one long trajectory. 
We start from a random particle configuration with uni- 
formly distributed orientations. After a sufficient large 
relaxation time (t — 25) we collect data for 500 time 
units. The coupling parameter F is then increased by 20 
and the protocol of relaxation and recording data is re- 
peated until we reach the maximal F. The melting runs 
for each pair {F, /} are independent starting out of the 
perfect hexagonal crystal albeit with random particle ori- 
entations. Again, we wait an adequate amount of time 
before collecting data for 50 time units. For both cooling 
and melting we record data from 5 independent runs for 
each{r,/}. 

We monitor structural changes through the global 
bond-orientational order parameter |25j 



1p6 



1 ^ 



(4) 



where Af{i) is the set of the six nearest neighbors of the 
ith particle and 9ij is the angle between the bond vector 
pointing from particle i to j and an arbitrary fixed axis. 
This order parameter is practically zero in the disordered 
phase whereas in a perfect crystal ipQ = 1. In Fig. [TJ^ we 
plot the -06 values averaged over all runs for both the cool- 
ing and the melting protocol. There is a clear transition 
between a disordered liquid and an ordered crystalline 
phase even for self-propelled particles (/ > 0). However, 
while the transition is rather abrupt for f — the struc- 
tural ordering is more gradual for higher propelling forces 
/. In Fig. we do not resolve a possible hexatic inter- 
mediate phase However, we note that no hysteresis 
is observed in agreement with a second-order transition 
scenario. As a structural criterion for both the melting 
and freezing transition we determine Fg from the condi- 
tion V'e = 0-45. In particular, for / = we find Fg ~ 240, 
which agrees well with a previous estimate [27) . 



Cooling the suspension, a dynamical criterion for freez- 
ing is given by the precipitous drop of the long-time dif- 
fusion coefficient 



D^ lim j-(|Ar,(t)p) 

t-i-oo 41 



(5) 



with Ar,{t) = r,(t) - ri{0). In Fig.[l^ we plot the diffu- 
sion coefficient for different forces. The value FJj at which 
the suspension freezes is estimated from the condition 
D = 0.086 ,23 . This gives an upper bound F < F£, to 
the liquid region, see the phase diagram Fig. [Tp. More- 
over, for not too large forces F^ ~ Fg correlates well 
with the position of the structural ordering as observed 
in Fig. Hence, this dynamical criterion for freezing 
based on particle mobility extends only to weakly driven 
suspensions of self-propelled particles. Note that at large 
forces and small F the diffusion coefficient D exceeds 1, 
the diffusion coefficient of a free passive Brownian parti- 
cle. 

We next consider a dynamical criterion for melting 
starting in the solid state and decreasing F. In one of 
the first theories for melting Lindemann conjectured that 
melting is caused by atom vibrations that start to in- 
terpenetrate [55]. In our case it is natural to consider 
the vibrational displacements of particles with respect 
to their lattice position. The Lindemann criterion then 
states that melting commences once the vibrational dis- 
placements reach a certain fraction of the lattice spac- 
ing. However, in two dimensions fluctuations on long 
wavelengths eventually destroy long-ranged positional or- 
der in the crystal |29j. The mean-squared displacement, 
therefore, is not a good measure to distinguish the liquid 
from the crystal. Instead, one defines a Lindemann-like 
parameter [3TJ [52] 



7l(0 



(|Ar.(t)-Ar,-(Op 
2P 



(6) 



from the neighbor-neighbor displacements. Here, i and j 
denote two particles that are initially neighbors. The lat- 
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FIG. 2: (Color online) (A) Time-dependence of the Linde- 
mann parameter Eq. Q for / = 8 below and at the melting 
point Fl — 700. (B) The plateau values of the Lindemann pa- 
rameter 7l both measured (closed symbols) and from Eq. ^ 
(open symbols) as a function of applied force. 



tice spadng of the hexagonal crystal is £ = 2^/^3~^/"' ~ 
1.075. In the liquid jL{t) diverges for long times with- 
out a plateau, whereas in the solid one observes a well 
defined plateau with Lindemann parameter 7l. Hence, 
we determine the melting point from the smallest value 
Fl for which we still observe a plateau with value 7l, 
see Fig. [2K. Above F > FJ^ the suspension is crystalline 
both with respect to orientational order and a vanishing 
diffusion coefficient. 

Using a simplified picture to describe the process of 
melting we assume particles to move independently close 
to their lattice position. The linearized forces then read 
—ViU « —k{Yi~Y^) with effective curvature A: oc F. The 
initial positions r,;(0) — r" correspond to lattice positions 
in the hexagonal crystal. A straightforward calculation 
of Eq. ([6| in the limit i — oo yields 



2p 



fc2 - £>? 



1 - 



k 



(7) 



In Fig. [2j3 the plateau value 7^ as a function of force 
is plotted together with the prediction 7l(/, rj^) from 
Eq. ([7]). Both values show excellent agreement. For the 
plot we have fixed the proportionality between k and F 
such that the values for f — are equal. Moreover, 
7£(0) ~ 0.026 agrees well with previous experiments [22]. 

While structural and dynamical criteria agree in equi- 
librium the phase diagram Fig. [Tp shows that in non- 
equilibrium there is a transition region Fq < F < F^ be- 
tween liquid and crystal, which widens for larger forces 
/. This region of parameter space is characterized by 
a high structural order but non-vanishing long-time dif- 
fusion. Moreover, the dynamical freezing and melting 
lines do not follow the orientational order but are shifted 
to higher F at higher forces. This implies that at high 
propelling speeds structural ordering occurs before dy- 
namical freezing. While an effective temperature could 
be defined individually for each criterion, the resulting 
values as a function of force clearly do not agree. 
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FIG. 3: (Color online) Top: Probability distributions for 
&t f — (solid lines) and f — 8 (dashed lines) for three 
different global values. Bottom: Snapshots of particle 
configurations for both the equilibrium (/ = 0, left column) 
and driven (/ = 8, right column) suspension. The rows cor- 
respond to constant global tp^ values: from top to bottom 
ip6 ~ 0, 0.45, 0.67, 0.8, cf. Fig. [Tp. Particles are colored ac- 
cording to their qe value. While liquid and crystal (top and 
bottom row) are indistinguishable the transition region (mid- 
dle rows) is marked by heterogeneous structure. 
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To gain further insight we define the order parameter 
geW^Re^ J2 ^mU) (8) 

per particle in order to effectively describe the local en- 
vironment of every particle [30 . The advantage of the 
neighbor-shell averaging compared to jggp is that more 
sharply distinguishes liquid-like from ordered regions. In 
Fig. [3] probability distributions of qe for / = and f — 8 
are plotted. For the crystal {tpe — 0.8) no difference 
between the driven and the undriven suspension is dis- 
cernible (see also the last row of Fig. |3]). In the liquid 
(tpQ ~ 0) the driven suspension is somewhat less struc- 
tured compared to equilibrium. This is caused by the 
larger effective diffusion due to the propulsion. A large 
difference can be seen in the distributions for suspensions 
with ipe ~ 0.45, i.e., in the transition regime. Here the 
driven suspension is locally more ordered but with a long 
tail that extends down to unordered particles. The spa- 
tial distribution of order and disorder corresponding to 
this ■06 for a single snapshot is shown in the second row of 
Fig. |3] For / = 8 the suspension is overall more ordered 
but also more heterogeneous, i.e., small, well separated 
liquid "bubbles" remain. Interestingly, the ipe = 0.67 
iso-line crosses the melting line such that for / = 8 it 
is within the transition region. Two snapshots for this 
case are depicted in the third row of Fig. [3j Due to the 
crossing the two forces now also describe two different 
dynamic regimes: while diffusion has effectively ceased 



in equilibrium, some particles still move in the driven 
suspension. 

In conclusion, we have shown by using Brownian dy- 
namics computer simulations that self-motile colloidal 
particles crystallize at sufficiently high densities. As 
compared to the equilibrium freezing of passive parti- 
cles there is a significant shift in the freezing density and 
additional large structural fluctuations appear caused by 
the self-propulsion. In principle, our predictions are veri- 
fiable in real-space experiments on colloidal model swim- 
mers on a (quasi) two-dimensional substrate [12[ llSj . 

In future work, it would be interesting to generalize 
our model to one which embodies an explicit swarming 
behavior such as a self-propelled rod model [5]. Fur- 
thermore, since equilibrium freezing is different in two 
and three spatial dimensions, it would be very interest- 
ing to simulate a corresponding three-dimensional model. 
Last but not least the influence of self-motility of the 
glass transition has not yet been studied. Since glass 
formation competes with crystallization and is typically 
accompanied with dynamical heterogeneity [5TI - I55] . self- 
propulsion may introduce an internal source of additional 
fluctuations which can help to form amorphous structures 
provided the density is large enough. 
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